Long term black hole evolution with the BSSN system by pseudospectral methods 
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We present long term evolutions of a single black hole of mass M with the BSSN system using 
pseudospectral methods. For our simulations we use the SGRID code where the BSSN system is 
implemented in its standard second order in space form. Previously we found that such simulations 
are quite unstable. The main goal of this paper is to present two improvements which now allow 
us to evolve for longer times. The first improvement is related to the boundary conditions at the 
excised black hole interior. We now use a gauge condition that ensures that all modes are going 
into the black hole, so that no boundary conditions are needed at the excision surface. The second 
more significant improvement has to do with our particular numerical method and involves filters 
based on projecting the double Fourier expansions used for the angular dependence onto Spherical 
Harmonics. With these two improvements it is now easily possible to evolve for several thousand 
M. The only remaining limitation seems to be the radiative outer boundary conditions used here. 
Yet this problem can be ameliorated by pushing out the location of the outer boundary, which leads 
to even longer run-times. 

PACS numbers: 04.25.dg, 04.30.Db, 04.70.Bw, 95.30.Sf 



I. INTRODUCTION 

Currently several gravitational wave detectors such as 
LIGO [l| or GEO600 0] are already operating, while sev- 
eral others are in the planning or construction phase [|[ . 
One of the most promising sources for these detectors are 
the inspirals and mergers of binary black holes. In order 
to make predictions about the final phase of such inspirals 
and mergers, fully non- linear numerical simulations of the 
Einstein Equations are required. In order to numerically 
evolve the Einstein equations, at least two ingredients 
are necessary. First we need a specific formulation of the 
evolution equations. And second, a particular numerical 
method such as finite differencing or a spectral method 
is needed to implement these equations on a computer. 
Currently there are two systems that have been used suc- 
cessfully to evolve binary black hole systems. There is 
the generalized harmonic system Q , which in its original 
second order in space form has been only used in a fi- 
nite differencing code 0, H, H, 0, IH • A nrs t order version 
of this system has also been successfull y ev olved using a 
spectral method 0, M, EH, E3, 03 Hill- The second 
and by far the most common system is the BSSN sys- 
tem This system is second order in space. Over the 
last few years, many successful binary black hole evolu- 
tions[l7Ul3Jl^ M, [30 , 
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However, long term simulations with 



with this system 
BSSN have so far only been performed with finite differ- 
encing techniques. It seems thus natural to ask whether 
the BSSN system in its standard second order form can 
also be implemented using a spectral method. This ques- 
tion has been partially answered in Ref. (5f| (henceforth 
paper 1) where we have found that a single black hole can 
in principle be evolved using BSSN together with a spec- 
tral method. However, we found that our simulations fail 
after a short time, usually after about lOOAf, where M is 



the mass of the black hole. This result has been obtained 
by using the same gauge conditions and outer boundary 
conditions for the spectral method as in successful finite 
difference implementations. In this paper we describe 
two new ingredients which allow us to evolve for much 
longer times with our spectral method. The first one is a 
gauge condition that ensures that no modes are leaving 
the black hole at the horizon, so that no boundary con- 
dition is needed at the horizon. The second ingredient 
is a particular spectral filter that is based on projecting 
the double Fourier expansions (used for the angular de- 
pendence of all fields) onto Spherical Harmonics. This 
filter is very useful in removing unphysical modes in our 
evolved fields and is the main reason for the observed 
increased run-times. 

Throughout we will use units where G = c = 1. The 
paper is organized as follows. In Sec. |TT] we describe 
briefly how a single black hole can be evolved with the 
SGRID code [H|,[56[ and we introduce the two new ingre- 
dients that lead to improved run-time. Sec. IIIII presents 
the results of several simulations using these ingredients. 
In Sec. IIVI we show results from two code tests in a 3- 
torus. We conclude with a discussion of our results in 
Sec. El 



II. BSSN AND THE SGRID CODE 

In this section we describe how a single black hole can 
be evolved using the SGRID code. We first discuss the 
BSSN system and then the spectral method we are using. 



A. A single black hole and the BSSN system 

As initial data we use a Schwarzschild black hole of 
mass M in Kerr-Schild coordinates. Thus initially the 
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3-metric and extrinsic curvature are given by 

'/,, = S ij + 2Hl i l j (1) 
Kij = a[kHj + I , II, + III,,, + Hl jti 

+2H 2 (lilklj,k + ljlkh,k) + 2HhljlkH,k], (2) 

where H = M/r and P = x l jr. The initial lapse and 
shift are 



a = 1/VT+2H, 
ft = 2l i H/{l + 2H). 



(3) 
(4) 



The BSSN formulation introduces a conformal met- 
ric 7ij (with determinant one) and conformal factor 
e 4 ^", so that = e^jij. It also introduce the addi- 
tional variable T l = j^j^jjkj- The extrinsic curva- 
ture is split into its trace K and tracefree part A^ using 
= e i(t> (Aij + #7y/ 3 )- The BSSN evolution equa- 
tions are then 



dtlij = 
d t <j> = 

d t r = 
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Here the superscript TF in Eq. ([8|) denotes the trace 
free part and Di is the derivative operator compatible 
with the 3-metric. In order to ensure that Aij remains 
traceless during our numerical evolution, we subtract any 
trace due to numerical errors after each evolution step 
from A^ . 

As in paper 1 we will keep the shift j3 l constant during 
all our evolutions and choose some particular coordinate 
condition to evolve the lapse a. Furthermore we employ 
spherical coordinates to evolve all fields inside a spheri- 
cal shell that is bounded by a inner radius Ri n located 
just inside the horizon and the outer radius R ou t that 
can be chosen freely. This means that we need bound- 
ary conditions at both the excision radius Ri n and at the 
outer boundary R ou t- At R ou t we use the same simple 
radiative boundary conditions that were also employed 
in paper 1, and which are used by virtually all other 
BSSN implementations. The difference in the current 
work lies at the excision boundary Ri n . In paper 1 we 
had simply assumed that all modes at the excision radius 
Rin are going into the hole (i.e. are leaving the numer- 
ical domain), which implies that no boundary condition 



should be imposed at Ri n . However, as already shown 
by Bona et al. [57], the gauge conditions of the "1+log" 
type used in paper 1 introduce faster than light gauge 
speeds so that one can expect gauge modes to leave the 
hole. By analyzing the characteristic modes of a first 
order in space extension of the BSSN system Beyer and 
Sarbach [581 ] have explicitly shown that the assumption 
that all modes are going into the hole is only true for par- 
ticular gauge conditions and that in fact certain gauge 
modes do leave the black hole for the gauges investigated 
in paper 1. This result still holds for the second order in 
space BSSN system used here. The reason is as follows: 
The extension in [58| consists of merely introducing the 
additional derivative variables d k = I2d k (f>, d k ij = dkjij, 
A k — (d k a)/a [78j and their evolution equations without 
adding any constraints. Hence the mode speeds of this 
first order system are the same as one would obtain from 
the formalism of Gundlach and Martin-Garcia [5^, [6(3] 
for second order systems, in which one would express the 
corresponding modes in terms of longitudinal derivatives 
of <fi, 7y and a. 

Thus not imposing any boundary conditions at the ex- 
cision radius while using the gauges in paper 1 will cause 
instabilities and is one of the reasons why the runs in 
paper 1 failed so early on. To avoid this problem we now 
use 



d t a = (i l dia 



a 2 K 



Cn 
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as the evolution equation for the lapse. Here the constant 

Co = a(t = 0) 2 K(t = 0) - P l (t = 0)d % a(t = 0) (11) 

is computed from the initial values of a, ft 1 and K and is 
chosen such that d t a = initially. Using the equations 
for the characteristic speeds given in (5S| it is clear that 
all modes (including the gauge modes) of the first order 
BSSN version are going into the black hole inside the 
horizon when the lapse evolves according to Eq. (| 10[) . 
Therefore it is now indeed justified to not impose any 
boundary conditions at Ri n , which is precisely what we 
will do here. This observation carries over to the second 
order in space version of the BSSN system we are using 
here. However, this new gauge condition alone does not 
lead to a big improvement in the run-time of our code. 
We also need the filter algorithm described in the next 
subsection. 



B. Numerical methods employed in the SGRID 
code 



For the time integration of all evolved fields we use a 
fourth order accurate Runge-Kutta time integrator. In 
order to compute spatial derivatives we apply essentially 
the same pseudospectral collocation method as in paper 
1 but this time we add a particular filter algorithm that 
leads to big increases in run-time. For a thorough discus- 
sion of collocation methods and other possible spectral 
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methods in a more general setting we refer the reader to 
a recent review article by Grandclement and Novak [6l[ . 

As in paper 1 we use standard spherical coordinates 
and use Chebyshev polynomials in the radial direction 
and Fourier expansions in both angles. The collocation 
points are then given by 



Rin Rout ( 
- COS 



2 \Nr-l 

9 3 = ir{2j + l)/N e 



Rin ~t~ Rout 



ipk = 2nk/N v , 



(12) 

(13) 
(14) 



where < i < N r - 1, < j < Ng - 1, < k < N v - 1 
and N r , Ng, N v are the number of collocation points in 
each direction. Notice that both angles run from to 
27T, which is necessary to ensure the periodicity in both 
angles needed for Fourier expansions. In paper 1 we have 
found that this double covering causes no problems when 
we evolve scalar fields. Furthermore, simply removing 
the double covering does not improve the run-time for a 
simulation with BSSN. 

In order to avoid problems near the coordinate sin- 
gularities at the poles we evolve the Cartesian compo- 
nents of all BSSN fields. The spatial components of the 
evolved fields thus are 0, 7^ , T l , Aij, and K. We can ex- 
pand these fields in terms of spherical harmonics. For a 
spherically symmetric hole the modes contained in these 
fields are / = for <fi and K, I = and / = 1 for T l , 
I = 0, I = 1 and I = 2 for 7^ and Ay. This shows 
that scalar, vector and tensor like fields contain a dif- 
ferent number of /-modes. However, our double Fourier 
expansions of both angles do not make this distinction, 
and in addition they even allow for fields that could not 
be represented by Spherical Harmonics (because of the 
double covering). This motivates the following spectral 
filter algorithm, which we have found to be crucial to 
extend the run-time of our simulations. From the repre- 
sentation of each field u by its values at the collocation 
points we compute expansion coefficients in terms of dis- 
crete Spherical Harmonics. I.e. for each radial grid point 
we compute [62| 



c(r-i)/ 



2B-1 



2B 



3=0 



k=0 



(15) 

where the P[ a (x) are Associated Legendre polynomials, 
their weights are 



W j 



2 . f(2j + l)n" 

B-l 



fc=0 



B = Ng/A and we choose Ng to be a multiple of 4. The 
sum over k in Eq. (|15[) is a simple Fourier transform, 
while the sum over j is computed using routines from the 
S2kit package [62, 63|. The maximum I in this expansion 



is Imax = B — 1 = Ng/A — 1, and as usual we set c(ri)i Tl 
if m > l m ax- Note that the sum over k in Eq. ([H 
only runs up to 2B — 1 = N$/2 — 1, i.e. the only field 
values taken into account come from points with < 9 < 
7r. From the coefficients c(rj)/ m we can then recompute 
u(ri,9j,(pk) via the inverse discrete Spherical Harmonic 
transform. 

In terms of a double Fourier representation an expan- 
sion up to l max corresponds to 2l max +l Fourier modes for 
both the 9- and ^-directions. I.e. after recomputing the 
u(ri,9j,tfk) via the inverse transform only 2l max + 1 = 
Ng/2 — 1 Fourier modes remain. Since originally we had 
Ng Fourier modes, this operation filters out about half 
of all Fourier modes for the 0-direction. Similarly it fil- 
ters out all modes above k — 2l max = Ng/2 — 2 for the 
(^-direction. For our simulations we choose N v = 3Ng/A. 
In that case the above projection onto Spherical Har- 
monics filters the upper third of the Fourier modes for 
the (^-direction. We use this procedure to filter all tensor 
like fields (i.e. 7^ and Ay). For vector like fields (T l and 
(3 l ) we additionally set c(ri)i maxm = so that we filter 
one more /-mode. For scalars like fields (cf>, K and a). 
we set c(ri)i maxm — c(ri); max _i m = so that we filter 
yet one more /-mode. Notice that this filter algorithm 
also removes the double covering. We apply it after each 
evolution step. It removes any unphysical modes (i.e. 
the ones which cannot be represented by an expansion 
of Spherical Harmonics). It also removes high frequency 
modes that have been contaminated by aliasing due to 
the non-linear terms in the BSSN system. Since the sub- 
mission of this paper we have also successfully used this 
filter in scalar field evolutions [64j . 

When we apply the filter algorithm for the angu- 
lar directions detailed above our code runs noticeably 
longer. Since the BSSN system is non-linear, aliasing 
also to occurs in the radial direction. We can further 
improve the run-time by applying the standard Orszag 
2/3 rule [6^, [6(| in the radial direction: We first com- 
pute the coefficients of the Chebyshev expansion in 
the radial direction for each field and then set ai — 
for all i > 2N r /3, i.e. we filter out the top 1/3 of the 
Chebyshev modes. From the so filtered we recompute 
the values of all fields at the collocation points. When 
we apply both the radial and especially the angular filter 
algorithms we get a greatly improved run-time. 



III. EVOLVING A BLACK HOLE WITH THE 
SGRID CODE 

In this section we show the results from several sim- 
ulations of a single black hole. The Kerr-Schild initial 
data together with the choice of lapse and shift explained 
above describe a static black hole, so that during evolu- 
tion all fields should remain constant. However due to 
numerical errors this is not exactly the case. In fact we 
observe that the norms of all fields grow in time until 
the code finally fails. This growth is exponential as is 
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FIG. 1: This plot shows the ADM mass Madm for simula- 
tions with different radial resolutions. In each case R ou t = 
480M and (Ng,N v ) = (16,12). The code performs much 
better for higher resolutions. 



FIG. 3: This plot shows the L -norm of the Hamiltonian con- 
straint violation for different angular resolutions. In each case 
Rout = 480M and iV r = 193. We see that the results do not 
depend much on the angular resolution. 
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FIG. 2: This plot shows the Z/ 2 -norm of the Hamiltonian con- 
straint violation for different radial resolutions. In each case 
R out = 480M and (iV e , 7V^) = (16, 12). We observe geometric 
convergence with increasing resolution. 



typical for the numerical instabilities in such problems. 
Nevertheless, the time scale for this growth is now much 
longer than in paper 1. 

In all the simulation described in this paper the time 
step in our Runge-Kutta time integrator is set to 

At = Ar mm /2, (17) 

where Ar m i n is the distance between the two closest grid 
points in the radial direction. For the first set of simu- 
lations we use a spherical shell that extends from Ri n = 
1.85M to Rout — 480M as our numerical domain. The 
angular resolution is chosen to be (Ng,N v ) — (16,12), 
which corresponds to l max — 3 after applying our filters. 
Figures [T] and show how the ADM mass Madm and 



the Hamiltonian constraint 

H = R + K 2 - K\K\ (18) 

evolve in time for different radial resolutions (given in 
terms of N r ). In all cases examined, our code fails even- 
tually. The time when it fails corresponds to the end of 
the lines in Fig. [2j However, the run-time in each case is 
significantly longer than in paper 1. Furthermore, we see 
from both figures that the code runs longer for higher res- 
olutions. For example for N r — 257 we get a run-time of 
3244M, which is more than 30 times longer than in paper 
1. As is evident from Fig. [2] our code exhibits the geo- 
metric convergence expected from a spectral code. From 
Fig. [1] we see that the Madm (which should in principle 
be always equal to M) remains closer to M for longer 
with increased resolution. 

The simulations discussed so far were all for the same 
angular resolution. However, as one can see from Fig. [31 
our results are largely independent of angular resolutions 
as one would expect for a spherically symmetric problem. 
For this reason we will revert to (Ng, N v ) — (16, 12) from 
now on. 

The reader may wonder why in fact our code fails 
at late times. The reason lies with the radiative outer 
boundary conditions we use when we evolve the BSSN 
system. Figure [4] depicts the normalized Hamiltonian 
constraint, which is obtained from the Hamiltonian con- 
straint divided by the magnitude of all terms that enter 
Eq. (|18[) . Compared to H this normalization magnifies 
errors near the outer boundary, since far away from the 
black hole the individual terms in Eq. (|18p are all small. 
Figure |4] shows that Hamiltonian constraint violations 
which arc at first (t = 500M ) concentrated near the outer 
boundary gradually enter the numerical domain, until the 
simulation fails shortly after t — 2500 M. In principle this 
problem should not come as a surprise. The standard ra- 
diative outer boundary conditions impose conditions on 
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Normalized Hamiltonian constraint (R =480M, N =193) 
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FIG. 4: This plot shows the normalized Hamiltonian along 
the radial direction (6 — 9ir/16, tp = 0) for different times. 
The resolution is (N r ,Ng,N v ) = (193,16,12) and the outer 
boundary is at Rout = 480Af. We see how constraint viola- 
tions enter through the outer boundary and in time contam- 
inate the entire numerical domain. 



FIG. 5: This plot shows the ADM mass for runs with different 
Rout- In each case we use (Ng,N v ) = (16, 12) and adjust N r 
such that the resolution near the black hole is the same for 
each run. The lines end when the run fails. We observe that 
the run-time increases with increasing R ou t- 



all evolved fields, which implies that we are imposing 
conditions on all ingoing and outgoing modes. However, 
we really should only impose conditions on the ingoing 
modes. In paper 1 we have shown that for a scalar field 
system the radiative conditions are equivalent to impos- 
ing conditions only on incoming modes, because the con- 
ditions on the outgoing modes are equivalent to the evo- 
lution equations. For the BSSN system, however, there is 
no clear cut answer. To our knowledge the well-posedness 
of the boundary problem for the BSSN system has never 
been established and remains an open question. The 
work by Beyer and Sarbach [58| mentioned above comes 
closest to addressing our problem. They provide explicit 
boundary conditions for the incoming modes of the BSSN 
system for the case of a frozen shift that is tangential to 
the boundary. These boundary conditions lead to a well 
posed initial-boundary value formulation. One problem 
is that our shift is not tangential to the outer bound- 
ary so that their proof does not strictly apply here. Yet 
even if the proof could be extended to a radial shift, an- 
other possibly more serious problem remains. Beyer and 
Sarbach provide only seven boundary conditions because 
there are exactly seven modes that enter through the 
outer boundary However, the radiative boundary 

conditions we use here impose conditions on all evolved 
variables, which amounts to 18 conditions. Simple count- 
ing arguments make it seem unlikely that our 18 condi- 
tions reduce to only 7 conditions on the incoming modes. 
So if indeed there is a problem with well-posedness at the 
continuum level, the real surprise is that radiative outer 
boundary conditions work as well as they do for the BSSN 
system. Recall that all the finite differencing codes that 
have recently been used to perform long term simulations 
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|4l|, |4J, |43|, |4J, |45|, |4fl |47j, |4& |49|, |5fl M m also use 
these radiative outer boundary conditions. These finite 
differencing codes are typically run for several thousand 
M without failing. However, they also have problems at 
the outer boundary in the sense that constraint viola- 
tions are observed to enter. So we conclude that we need 
improved outer boundary conditions if we want to evolve 
for even longer times, at least with a spectral code. Nev- 
ertheless, an evolution time of about 2500M is not all 
that bad and may already be sufficient for many prob- 
lems. In addition, we have found a simple way to further 
increase the run-time. Figure [5] shows the ADM mass 
Madm for runs with different R ou t- The lines end when 
the code fails. In all these runs we use the same angular 
resolution and we adjust N r such that the radial resolu- 
tion is kept constant near the black hole horizon. We see 
that the run-time is longer for larger R ou t- Thus we are 
able to increase the run-time as much as we need. The 
only cost is the additional computer time needed to sim- 
ulate a larger domain. For example for R ou t = 1920A/ 
we can evolve for more than 4500M. Note that while the 
radial resolution is the same near the black hole horizon 
for all runs shown in Fig. [5l the resolution in the mid- 
dle of the domain is less for the runs with larger R ou t- 
This fact is related to the uneven spread of collocation 
points for the Chebyshev expansions used in the radial 
direction. Thus effectively the runs with larger R ou t are 
less accurate. We could have of course increased N r even 
further so that the resolution in the middle of the domain 
would always be the same. In that case the radial resolu- 
tion near the black hole would increase with R ou t and we 
would get even longer run-times, since as evidenced by 
Fig. [2] the run-time only increases with increased resolu- 
tion. Figure [HI is the analog of Fig. [5] for the case where 
we increase N r proportional to R ou t- We see that the 
run-time is now approximately proportional to R ou t ■ 
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FIG. 6: This plot shows the Hamiltonian constraint for runs 
with different Rout- In each case we use (N$,N V ) = (16, 12) 
and adjust iV r such that the resolution in the middle of the 
domain is the same for each run. The lines end when the run 
fails. We observe that the run-time is proportional to R ou t- 
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FIG. 8: The normalized Hamiltonian in the xz-pla.ne at 
t = 40M for a run using our spectral filters. The resolution 
is (N r , No, Nip) — (33,16,12) and the outer boundary is at 
Rout = 60M. We see how constraint violations have entered 
through the outer boundary. 
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FIG. 7: This plot shows the lapse a along the radial di- 
rection for the run with Rout = 480AZ and (N r , N e , N v ) = FIG. 9: Same as Fig.® but at t = 88M. 
(257, 16, 12) for three different times. 



In order to get a sense of how much the metric fields 
change during our evolutions Fig. [7] depicts the lapse a at 
three different times for the longest run of Fig. [6] From 
Fig. [7] we see that until t = 3000M the lapse does not 
change by much. Most of the change occurs just before 
the run fails at 3176M. 

In order to study the angular dependence of the con- 
straint violations entering through the outer boundary 
we next show snapshots of the normalized Hamiltonian 
constraint in the ccz-plane. For visual clarity we 

depict results from a short low resolution run with outer 
boundary at R out — 60M. Figures ISlfTTl show that the 
constraint violations have no strong angular dependence 
when our filters are active. Yet we see strong oscillations 
in the radial direction. This particular simulation fails 
shortly after t = 424M. Figures IT21 and [T3l show how our 
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FIG. 10: Same as Fig.® but at t = 200A/. 
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Hamiltonian constraint 




FIG. 11: Same as Fig. E] but at t = 424M just before the run 
fails. 




results change if our filters are switched off. The con- 
straint violations are larger, the run fails sooner (shortly 
after t = 88M), and towards the end our results deviate 
significantly from spherical symmetry. Nevertheless the 
strongest oscillations are still in the radial direction. 



IV. BSSN TESTS IN A 3-TORUS 

As we have seen our black hole runs always fail at 
some point. In order to strengthen the conjecture that 
this failure is related to the radiative outer boundary 
conditions, we have also performed two simple tests in a 
3-torus, i.e. without boundary conditions. These tests 
were originally proposed in b 67] . For the BSSN system 
both tests have be en p erformed before [68| using a finite 
difference code |69j, llOj • As we will see below the spectral 
SGRID code is at least as stable as the finite differencing 
code. Both tests are carried out in Cartesian coordinates 
using Fourier expansions in x, y and z-directions. Also, 
both are effectively one-dimensional as we only use either 
3 or 1 points in the y and z-directions. 



A. Random perturbations on flat space 

In this test we evolve flat space in a 3-torus with added 
random perturbations. The initial data are given by [67j 

9ij = % + £ ij, (19) 

Ku = eij, (20) 

a = 1 + Eij, (21) 

P l = 0, (22) 



FIG. 12: Same as Fig. [8]but without filtering. 



Hamiltonian constraint 




FIG. 13: Same as Fig.[TJ]but at t = 88M just before the run 
fails. No niters were used. The result is no longer spherically 
symmetric. 



where ey is a random number with a proba- 
bility distribution that is uniform in the interval 
[-IQ- 10 /(50/N X ) 2 ,+1Q- 10 /(50/N X )% Hence our initial 
data differ slightly from the ones in [681 ] who add random 
noise only to gtj. The parameters for this test are: 

• Points: N x is varied, N y = N z = 3 

• Resolution Ax = 1/N X 

• Simulation domain: x G [0, 1], (y, z) € [0, 3Ax] 

• Gauge: d t a = —a 2 trK, (3 l = 

Our results are shown in Fig. [14] As we can see our 
simulations are all stable. The Hamiltonian constraint 
violation is constant. As expected the Hamiltonian con- 
straint does not converge to zero, since we have added 
random perturbations to flat space, which do not sat- 
isfy the constraints. The amplitude of the added random 
noise was chosen such that the Hamiltonian constraint vi- 
olation should be independent of the spatial resolution. 
From Fig. [TJ]we see that we indeed get roughly the same 
constraint violation for each run. The same qualitative 
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FIG. 14: The L -norm of the Hamiltonian constraint does 
not grow in time if we evolve fiat space with random noise 
in a 3-torus. As expected, the Hamiltonian constraint has 
approximately the same value for all resolutions. 



FIG. 15: The L -norm of the Hamiltonian constraint grows 
exponentially for a gauge wave in a 3-torus. Convergence is 
lost after about 600 crossing times. The lines end when the 
runs fail. Nevertheless the runs last about seven times longer 
than with a finite difference code. 



results for BSSN were obtained in [68| with a finite dif- 
ferencing code. We should also note that this test is non- 
trivial. In 68] it was shown that (at least in the finite 



differencing case) both the ADM [71| and the AA [72 
system fail this test. 



B. Gauge wave 



method. Furthermore, the gauge wave test is a member 
of a one-parameter family of exponential harmonic gauge 
solutions [73l . l74| . Any numerical error can probe a mem- 
ber that is exponentially growing. Thus generically we 
expect an exponential blow up for any formulation with 
any code. 



In this test we look at the stability of the BSSN sys- 
tem in the non-linear case. We consider flat space in a 
coordinate system where the metric is given by [67] ] 

ds 2 = -(1 + a)dt 2 + (1 + a)dx 2 + dy 2 + dz 2 , (23) 

with a — A sin [2 (x — t)/L]. We set our initial data using 
this metric at t = 0. Here L is the size of the domain in 
the ir-direction and A is the amplitude of the wave. Since 
this wave propagates along the a;-axis and all derivatives 
are zero in the y- and z-directions, the problem is one- 
dimensional. For our tests we use these parameters: 

• Points: N x is varied, N y = N z = 1 

• One-dimensional simulation domain: x € [0, 1] 

• A = 0.01 

• Gauge: d t a = —a 2 trK, /?* = 

From Fig. [TS] we see that our simulations are unstable. 
We lose convergence after about 600 light crossing times. 
After that the runs fail at around 700 crossing times. 
While this result is disappointing, it is not surprising. 
When Jansen et al. [68[ performed the very same test 
with a finite differencing code, BSSN failed even earlier 
at around 100 crossing times. So in fact the SGRID code 
runs longer than a finite differencing code and the ob- 
served instability is likely not a result of our numerical 



C. Test results 

We should point out that for the two test results pre- 
sented above we have used the Orszag 2/3 rule in the x- 
direction to filter out high frequency Fourier modes. The 
motivation for this filter was the success of the filters in 
our black hole evolution. We have found, however, that 
this filter does not help and that we obtain qualitatively 
the same results without it. 

As we have seen, the spectral implementation of BSSN 
performs as well or better than a finite differencing im- 
plementation. In the random noise test the BSSN imple- 
mentation in SGRID shows long term stability. In the 
gauge wave test SGRID runs longer than a finite differ- 
encing code, but fails after about 700 crossing times. We 
note that our black hole simulations last only a few cross- 
ing times. Thus when measured in crossing times BSSN 
in SGRID runs much longer without boundary conditions 
than with radiative boundary conditions. We interpret 
this fact as further evidence for our conjecture that the 
radiative outer boundary conditions are the main reason 
for the failure of our black hole evolutions. 



V. DISCUSSION 

As we have seen long term evolutions of single black 
holes with the BSSN system are possible with a spec- 
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tral method. One of the ingredients needed to achieve 
this goal is a suitable gauge condition that is compati- 
ble with imposing no boundary conditions at the excision 
surface. The second and main ingredient is the spectral 
filter described in Sec. Ill Bl above. This filter is capa- 
ble of removing unphysical modes related to the double 
covering introduced by double Fourier expansions. At 
the same time it also removes high frequency modes that 
are contaminated by aliasing due to non-linear terms in 
the BSSN system. Simulations with this filter run sig- 
nificantly longer. The only limitation seem to be prob- 
lems at the outer boundary coming from the simplistic 
radiative conditions used there. This shows that further 
work along the lines of [75[ is in principle needed to de- 
termine better outer boundary conditions for the BSSN 
system. There has been some work [5^, H3, [ZSl on de- 
termining ingoing and outgoing modes for the BSSN sys- 
tem in its standard second order in space form. Also, 
a recent paper by Nunez and Sarbach [77| presents in- 
teresting new boundary conditions that lead to a well 
posed initial-boundary value formulation for the case of 
linearized gravity. However, we are not aware of any pub- 
lication that provides explicit boundary conditions for 
BSSN which we could directly implement in our code. 
And deriving such boundary conditions is beyond the 
scope of this work. But as we have seen, the problem 
of not having good outer boundary conditions can also 
be circumvented by pushing out the location of the outer 
boundary. In addition, a run-time of several thousand M 
may be completely sufficient for many applications. For 
example, the many simulations we recently performed 
with our finite differencing code to determine the final 
mass and spin of black hole mergers [47| all just needed 



a run-time between 400M and 1000M. So it seems that 
the main roadblock for long term simulations of the BSSN 
system with spectral methods has been removed by the 
filter algorithms introduced in Sec. Ill Bl 

Nevertheless, our work provides a strong motivation for 
finding better outer boundary conditions for the BSSN 
system. We have observed that radiative outer boundary 
conditions cause constraint violations to enter the numer- 
ical domain. This happens for both spectral and finite 
differencing codes. Furthermore, BSSN in our spectral 
code runs for only a few light crossing times when we 
evolve with radiative outer boundary conditions, while it 
runs for hundreds of crossing times without these bound- 
ary conditions. So far this instability has not been ob- 
served with finite differencing codes. One reason for this 
difference may be that finite differencing codes are lo- 
cal since their stencils cover only a few points. Spectral 
codes, on the other hand, are global in the sense that 
their stencils cover all points, so that problems at the 
boundary will immediately affect all points. Thus they 
are inherently more sensitive to problematic boundary 
conditions. However, it is possible that the boundary 
driven instabilities we observe with SGRID could also af- 
fect finite differencing codes in some situations, possibly 
at later times. 

Acknowledgments 

It is a pleasure to thank Peter Diener and David 
Hilditch for useful discussions about BSSN gauge modes 
near black holes and boundary conditions. This work 
was supported by NSF grants PHY-0652874 and PHY- 
0855315. 



[1] LIGO, http://www.ligo.caltech.edu/. 
[2] GEO600, http://www.geo600.org/. 
[3] B. Schutz, Class. Quantum Grav. 16, A131 (1999). 
[4] F. Preterms, Class. Quant. Grav. 22, 425 (2005), gr- 
qc/0407110. 

[5] F. Preterms, Phys. Rev. Lett. 95, 121101 (2005), gr- 
qc/0507014. 

[6] F. Preterms, Class. Quant. Grav. 23, S529 (2006), gr- 
qc/0602115. 

[7] F. Pretorius and D. Khurana, Class. Quant. Grav. 24, 

S83 (2007), gr-qc/0702084. 
[8] C. Palenzuela, M. Anderson, L. Lehner, S. L. Liebling, 

and D. Neilsen, Phys. Rev. Lett. 103, 081101 (2009), 

0905.1121. 

[9] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, 
and O. Rinne, Class. Quant. Grav. 23, S447 (2006), gr- 
qc/0512093. 

[10] M. A. Scheel et al., Phys. Rev. D74, 104006 (2006), gr- 
qc/0607056. 

[11] M. Boyle et al, Phys. Rev. D76, 124038 (2007), 
0710.0158. 

[12] M. A. Scheel et al., Phys. Rev. D79, 024003 (2009), 
0810.1767. 



[13] E. Pazos, M. Tiglio, M. D. Duez, L. E. Kidder, and S. A. 

Teukolsky, Phys. Rev. D80, 024027 (2009), 0904.0493. 
[14] B. Szilagyi, L. Lindblom, and M. A. Scheel (2009), 

0909.3557. 

[15] T. Chu, H. P. Pfeiffer, and M. A. Scheel (2009), 
0909.1313. 

[16] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, 

024007 (1998), gr-qc/9810065. 
[17] M. Campanelli, C. O. Lousto, P. Marronetti, and 

Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr- 

qc/0511048. 

[18] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr- 
qc/0511103. 

[19] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. D73, 104002 (2006), gr- 
qc/0602026. 

[20] J. G. Baker et al., Astrophys. J. 653, L93 (2006), astro- 
ph/0603204. 

[21] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. 

Rev. D73, 061501 (2006), gr-qc/0601091. 
[22] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. 

Rev. D74, 041501 (2006), gr-qc/0604012. 



10 



[23] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. 

Rev. D74, 084023 (2006), astro-ph/0608275. 
[24] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krish- 

nan, and D. Merritt, Phys. Rev. D75, 064030 (2007), 

gr-qc/0612076. 

[25] J. A. Gonzalez, U. Sperhake, B. Bruegmann, M. Han- 
nam, and S. Husa, Phys. Rev. Lett. 98, 091101 (2007), 
gr-qc/0610154. 

[26] U. Sperhake, Phys. Rev. D76, 104015 (2007), gr- 
qc/0606079. 

[27] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer- 
ritt, Phys. Rev. Lett. 98, 231102 (2007), gr-qc/0702133. 

[28] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. 
Rev. D77, 101501 (2008), 0710.0879. 

[29] J. A. Gonzalez, M. D. Hannam, U. Sperhake, B. Brueg- 
mann, and S. Husa, Phys. Rev. Lett. 98, 231101 (2007), 
gr-qc/0702052. 

[30] B. Bruegmann, J. A. Gonzalez, M. Hannam, S. Husa, and 
U. Sperhake, Phys. Rev. D77, 124047 (2008), 0707.0135. 

[31] F. Herrmann, I. Hinder, D. M. Shoemaker, P. Laguna, 
and R. A. Matzner, Phys. Rev. D76, 084032 (2007), 
0706.2541. 

[32] I. Hinder, B. Vaishnav, F. Herrmann, D. Shoemaker, and 
P. Laguna, Phys. Rev. D77, 081502 (2008), 0710.5167. 

[33] M. Koppitz et al., Phys. Rev. Lett. 99, 041102 (2007), 
gr-qc/0701163. 

[34] P. Marronetti, W. Tichy, B. Bruegmann, J. Gonzalez, 

M. Hannam, S. Husa, and U. Sperhake, Class. Quant. 

Grav. 24, S43 (2007), gr-qc/0701123. 
[35] P. Marronetti, W. Tichy, B. Brugmann, J. Gonzalez, and 

U. Sperhake, Phys. Rev. D77, 064010 (2008), 0709.2160. 
[36] D. Pollney et al., Phys. Rev. D76, 124002 (2007), 

0707.2559. 

[37] L. Rezzolla et al., Astrophys. J679, 1422 (2008), 
0708.3999. 

[38] L. Rezzolla et al., Astrophys. J. 674, L29 (2008), 
0710.3345. 

[39] L. Rezzolla et al., Phys. Rev. D78, 044002 (2008), 
0712.3541. 

[40] U. Sperhake et al., Phys. Rev. D78, 064069 (2008), 
0710.3823. 

[41] W. Tichy and P. Marronetti, Phys. Rev. D76, 061502 

(2007) , gr-qc/0703075. 

[42] S. Dain, C. O. Lousto, and Y. Zlochower, Phys. Rev. 

D78, 024039 (2008), 0803.0351. 
[43] B. Brugmann et al., Phys. Rev. D77, 024027 (2008), gr- 

qc/0610128. 

[44] J. Healy et al., Phys. Rev. Lett. 102, 041101 (2009), 
0807.3292. 

[45] I. Hinder, F. Herrmann, P. Laguna, and D. Shoemaker 

(2008) , 0806.1037. 

[46] C. O. Lousto and Y. Zlochower, Phys. Rev. D79, 064018 

(2009) , 0805.0159. 

[47] W. Tichy and P. Marronetti, Phys. Rev. D78, 081501 

(2008), 0807.2985. 
[48] M. C. Washik et al., Phys. Rev. Lett. 101, 061102 (2008), 

0802.2520. 

[49] T. Bode et al., Phys. Rev. D80, 024008 (2009), 
0902.1127. 

[50] H. Nakano, M. Campanelli, C. O. Lousto, and Y. Zlo- 
chower (2009), 0901.3861. 

[51] B. Aylott et al., Class. Quant. Grav. 26, 165008 (2009), 
0901.4399. 



[52] B. Aylott et al., Class. Quant. Grav. 26, 114008 (2009), 
0905.4227. 

[53] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, and 

P. Diener (2009), 0910.3803. 
[54] C. Reisswig et al. (2009), 0907.0462. 
[55] W. Tichy, Phys. Rev. D74, 084005 (2006), gr- 

qc/0609087. 

[56] W. Tichy, Class. Quant. Grav. 26, 175018 (2009), 
0908.0620. 

[57] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 

Lett. 75, 600 (1995), gr-qc/9412071. 
[58] H. Beyer and O. Sarbach, Phys. Rev. D70, 104004 

(2004), gr-qc/0406003. 
[59] C. Gundlach and J. M. Martin-Garcia, Phys. Rev. D70, 

044031 (2004), gr-qc/0402079. 

[60] C. Gundlach and J. M. Martin-Garcia, Phys. Rev. D70, 

044032 (2004), gr-qc/0403019. 

[61] J. N. Philippe Grandclement, Living Re- 
views in Relativity 12 (2009), URL 
http : //www . livingreviews . org/lrr- 2009-1 

[62] P. K. D. Healy Jr., D. Rockmore and S. Moore, Journal 
of Fourier Analysis and Applications 9(4), 341 (2003). 

[63] P. J. Kostelec and D.N. Rockmore, S2kit: A Lite Version 
of SpharmonicKit: 

http : //www . cs .dartmouth.edu/ ~geelong/ sphere/. 

[64] I. Vega, P. Diener, W. Tichy, and S. Detweiler, Phys. 
Rev. D80, 084021 (2009), 0908.2138. 

[65] S. A. Orszag, J. Atmosph. Sci. 28, 1074 (1971). 

[66] J. P. Boyd, Chebyshev and Fourier Spectral Methods (Sec- 
ond Edition, Revised) (Dover Publications, New York, 
2000). 

[67] M. Alcubierre, G. Allen, T. W. Baumgarte, C. Bona, 
D. Fiske, T. Goodale, F. S. Guzman, I. Hawke, S. Hawley, 
S. Husa, et al., Class. Quantum Grav. 21, 589 (2004), gr- 
qc/0305023. 

[68] N. Jansen, B. Bruegmann, and W. Tichy, Phys. Rev. 
D74, 084022 (2006), gr-qc/0310100. 

[69] B. Bruegmann, W. Tichy, and N. Jansen, Phys. Rev. 
Lett. 92, 211101 (2004), gr-qc/0312112. 

[70] B. Bruegmann, J. Gonzalez, M. Hannam, S. Husa, 
U. Sperhake, and W. Tichy, Phys. Rev. D77, 024027 
(2008), gr-qc/0610128. 

[71] R. Arnowitt, S. Deser, and C. W. Misner, in Gravita- 
tion: An Introduction to Current Research, edited by 
L. Witten (John Wiley, New York, 1962), pp. 227-265, 
gr-qc/0405109. 

[72] A. Alekseenko and D. Arnold, Phys. Rev. D68, 064013 

(2003), gr-qc/02 10071. 
[73] M. C. Babiuc, B. Szilagyi, and J. Winicour, Class. Quant. 

Grav. 23, S319 (2006), gr-qc/0511154. 
[74] V. Paschalidis, J. Hansen, and A. Khokhlov, Phys. Rev. 

D78, 064048 (2008), 0712.1258. 
[75] G. Calabrese, J. Pullin, O. Sarbach, M. Tiglio, and 

O. Reula, Commun. Math. Phys. 240, 377 (2003), gr- 

qc/0209017. 

[76] C. Gundlach and J. M. Martin-Garcia, Class. Quant. 

Grav. 23, S387 (2006), gr-qc/0506037. 
[77] D. Nunez and O. Sarbach (2009), 0910.5763. 
[78] Any linear combination of these additional variables 

would lead to the same mode speeds. 
[79] This fact also holds for a shift that is not tangential. 



